function [price_coeff, elas_table, mc, w, brand_id_by_mkt, price_pre, share_pre, ...
    mean_mkup, mean_mkup_merge, price_cf, w_cf, ....
    mean_price_pre, mean_price_cf, mean_price_merge_pre, mean_price_merge_cf, ...
    mkt_weight, diver_craft, diver_outside, diver_lager, diver_ale]...
    = elas_mc_beer(para_nonlinear, ...
    para_linear, meanutility, price_data, mktid, x_rand, dim_price, ...
    randomdraw0, rand_income, income_norm, ...
    ind_mkt_start, ind_mkt_end, ...
    setup, mkt_set0, merge_craft, x_rand2, brand_id, craftdummy, ...
    oem_ind, mktsize)
% compute elasticities and invert the marginal cost


n_obs = size(price_data, 1);
para_nonlinear = para_nonlinear';

n_x_rand = size(x_rand, 2);
[randomdraw, para_income] = ...
    gen_randomdraw(para_nonlinear, n_x_rand, randomdraw0, false);

n_x_rand2 = size(x_rand2, 2);    
mu = [x_rand, ones(n_obs, 1)] * randomdraw;    
mu = mu + (para_income*rand_income).*price_data/income_norm ...
    + para_nonlinear(n_x_rand+4)*rand_income/income_norm ...
    + (x_rand2*para_nonlinear(n_x_rand+5:n_x_rand+n_x_rand2+4, :)).*rand_income/income_norm;
expmu = exp(mu); 

[share, prob_i] = ind_shnorm(exp(meanutility),expmu,ind_mkt_start,ind_mkt_end,mktid);

price_coeff = para_income*rand_income/income_norm + para_linear(dim_price);    
has_positive = true;    
while has_positive
    if any(price_coeff(:)>0)
        para_income = para_income*0.99;
        price_coeff = para_income*rand_income/income_norm + para_linear(dim_price);
    else
        has_positive = false;
    end
end    

% if merge_craft is not empty, redefine oem_ind to reflec the merge
if ~isempty(merge_craft)    
    oem_ind_cf = oem_ind>0;    
    oem_ind_cf(:, merge_craft(1)) = sum(oem_ind(:, merge_craft), 2)>0;    
    oem_ind_cf(:, merge_craft(2:end)) = 0;
end

n_mkt = max(mktid);

if isempty(mkt_set0)    
    mkt_set = 1 : n_mkt;    
else    
    mkt_set = reshape(mkt_set0, [1 length(mkt_set0)]);    
end

elas_table = cell(n_mkt, 1);
w = cell(n_mkt, 1);
share_pre = cell(n_mkt, 1);
mc = cell(n_mkt, 1);
price_pre = cell(n_mkt, 1);
price_cf = cell(n_mkt, 1);
w_cf = cell(n_mkt, 1);

mean_price_pre = nan(n_mkt, 1);
mean_price_cf = nan(n_mkt, 1);
mean_price_merge_pre = nan(n_mkt, 1);
mean_price_merge_cf = nan(n_mkt, 1);
mean_mkup = nan(n_mkt, 1);
mean_mkup_merge = nan(n_mkt, 1);

brand_id_by_mkt = cell(n_mkt, 1);

diver_craft = nan(n_obs, 1);% diversion ratios
diver_ale = nan(n_obs, 1);
diver_lager = nan(n_obs, 1);
diver_outside = nan(n_obs, 1);
mkt_weight = nan(n_mkt, 1);


carrier_ind = []; % do not consider retailer pricing

for nm = mkt_set    
    mkt_weight(nm) = mean(mktsize(mktid==nm));    
    prob_i_nm = prob_i(mktid==nm, :);    
    price_coeff_nm = price_coeff(mktid==nm, :);    
    n_nm = sum(mktid==nm);    
    prob_i_alpha_nm = prob_i_nm.*price_coeff_nm;    
    dsdp_nm = 0;
    
    for im = 1 : size(prob_i, 2)
        dsdp_nm = dsdp_nm - bsxfun(@times, prob_i_alpha_nm(:, im), prob_i_nm(:, im)')/size(prob_i, 2);
    end
    
    dsdp_nm = dsdp_nm + bsxfun(@times, eye(n_nm), mean(prob_i_alpha_nm, 2));
    
    craft_nm = craftdummy(mktid==nm);    
    lager_nm = x_rand(mktid==nm, 1);    
    ale_nm = x_rand(mktid==nm, end);   

    diver_craft(mktid==nm) = sum(dsdp_nm.*craft_nm.*(eye(n_nm)==0))'./abs(diag(dsdp_nm));
    diver_outside(mktid==nm) = 1 - sum(dsdp_nm.*(eye(n_nm)==0))'./abs(diag(dsdp_nm));    
    diver_lager(mktid==nm) = sum(dsdp_nm.*lager_nm.*(eye(n_nm)==0))'./abs(diag(dsdp_nm));    
    diver_ale(mktid==nm) = sum(dsdp_nm.*ale_nm.*(eye(n_nm)==0))'./abs(diag(dsdp_nm));
    elas_table{nm} = dsdp_nm.*(price_data(mktid==nm)'.*(1./mean(prob_i_nm, 2)));
    
    [markup_owner_nm_all, negativemarkup_retail_nm_all] = invertcost(share,prob_i,...
        price_coeff,...
        ind_mkt_end,ind_mkt_start, carrier_ind, oem_ind,...
        nm, setup);
    
    if isequal(setup.pricing_model, 'direct_oem')        
        mc{nm} = price_data(mktid==nm) - markup_owner_nm_all(mktid==nm);        
    elseif isequal(setup.pricing_model, 'sequential')        
        w{nm} = price_data(mktid==nm) + negativemarkup_retail_nm_all(mktid==nm);        
        mc{nm} = w{nm} - markup_owner_nm_all(mktid==nm);        
    else
        error('to be written');        
    end
    
    mkup_nm = price_data(mktid==nm & craftdummy) - mc{nm}(craftdummy(mktid==nm));    
    mean_mkup(nm) = sum(mkup_nm.*share(mktid==nm & craftdummy))/sum(share(mktid==nm & craftdummy));
    
    if ~isempty(merge_craft)
        mean_mkup_merge(nm) = sum(mkup_nm(oem_ind_cf(mktid==nm & craftdummy, merge_craft(1))).*share(mktid==nm & craftdummy & oem_ind_cf(:, merge_craft(1))))/sum(share(mktid==nm & craftdummy & oem_ind_cf(:, merge_craft(1))));
    end
    
    price_pre{nm} = price_data(mktid==nm);    
    share_pre{nm} = share(mktid==nm);
    
    if ~isempty(merge_craft)        
        n_prod_nm = sum(mktid==nm);
        same_firm_nm = false(n_prod_nm, n_prod_nm);
        same_firm_nm_cf = false(n_prod_nm, n_prod_nm);
        
        for k = 1 : size(oem_ind_cf,2)            
            ind_nm_k = oem_ind(mktid==nm, k)>0;
            
            if any(ind_nm_k)                
                same_firm_nm(ind_nm_k, ind_nm_k) = true;                
            end            
        end
        
        for k = 1 : size(oem_ind_cf,2)            
            ind_nm_k = oem_ind_cf(mktid==nm, k)>0;
            
            if any(ind_nm_k)                
                same_firm_nm_cf(ind_nm_k, ind_nm_k) = true;                
            end            
        end
        
        mean_wo_price = meanutility - para_linear(dim_price)*price_data;
        
        if isequal(setup.pricing_model, 'direct_oem')            
            setup.mode = 2;            
            [p_nm_test, ~, share_test] = equi_pwq(...
                para_nonlinear, randomdraw0, ...
                price_data(mktid==nm), price_data(mktid==nm)*0.8,...
                x_rand(mktid==nm, :), mean_wo_price(mktid==nm), true(n_prod_nm, 1), ... craftdummy(mktid==nm),...
                [], oem_ind(mktid==nm, :)>0, ...
                price_coeff(mktid==nm, :), rand_income(mktid==nm, :),...
                income_norm, mc{nm}, setup, x_rand2(mktid==nm, :), same_firm_nm);
            
            if max(abs(p_nm_test - price_pre{nm}))>1e-4
                error('mc inconsistent with direct oem pricing');
            end
            
            [p_nm_cf, ~, share_cf] = equi_pwq(...
                para_nonlinear, randomdraw0, ...
                price_data(mktid==nm), price_data(mktid==nm)*0.8,...
                x_rand(mktid==nm, :), mean_wo_price(mktid==nm), true(n_prod_nm, 1), ..., craftdummy(mktid==nm),...
                [], oem_ind_cf(mktid==nm, :)>0, ...
                price_coeff(mktid==nm, :), rand_income(mktid==nm, :),...
                income_norm, mc{nm}, setup, x_rand2(mktid==nm, :), same_firm_nm_cf);
            
            price_cf{nm} = p_nm_cf;       
            mean_price_pre(nm) = sum(p_nm_test(craftdummy(mktid==nm)).*share_test(craftdummy(mktid==nm)))/sum(share_test(craftdummy(mktid==nm)));            
            mean_price_cf(nm) = sum(p_nm_cf(craftdummy(mktid==nm)).*share_cf(craftdummy(mktid==nm)))/sum(share_cf(craftdummy(mktid==nm)));            
            ind_merge = craftdummy(mktid==nm) & oem_ind_cf(mktid==nm, merge_craft(1));            
            mean_price_merge_pre(nm) = sum(p_nm_test(ind_merge).*share_test(ind_merge))/sum(share_test(ind_merge));            
            mean_price_merge_cf(nm) = sum(p_nm_cf(ind_merge).*share_cf(ind_merge))/sum(share_cf(ind_merge));
            
        elseif isequal(setup.pricing_model, 'sequential')            
            setup.model = 1;            
            [p_nm_cf, w_nm_cf] = equi_pwq(...
                para_nonlinear, randomdraw0, ...
                price_data(mktid==nm), price_data(mktid==nm)*0.8,...
                x_rand(mktid==nm, :), meanutility(mktid==nm), craftdummy(mktid==nm),...
                carrier_ind(mktid==nm, :), oem_ind(mktid==nm, :), ...
                price_coeff(mktid==nm, :), rand_income(mktid==nm, :),...
                income_norm, mc{nm}, setup, x_rand2(mktid==nm, :), same_firm_nm);            
            price_cf{nm} = p_nm_cf;
            w_cf{nm} = w_nm_cf;
        else
            error('to be written');
        end
    end
    
    if ~isempty(brand_id)
        brand_id_by_mkt{nm} = brand_id(mktid==nm, :);
    end
end




function [meanprob, prob_i, exp_u] = ind_shnorm(expmeanval,expmu,ind_mkt_start,ind_mkt_end,mktid)
%% IND_SHNORM
% This function computes the "individual" probabilities of choosing each brand.
% The probabilities are those associated with the normally-distributed r.c. logit.
%
% ARGUMENTS:
% expmeanval = vector of exponentiated mean utilities, sorted by market then product
% expmu = matrix of exponentiated draws of deviations from mean utilities, sorted by market then product
%
% OUTPUT:
% meanprob = vector of expected market shares, sorted by market then product


exp_u = bsxfun(@times, expmeanval, expmu);

n_mkt = length(ind_mkt_start);

denom_exp_u = nan(n_mkt, size(expmu,2));
for t = 1 : n_mkt
    ind = ind_mkt_start(t):ind_mkt_end(t);
    denom_exp_u(t, :) = sum(exp_u(ind, :), 1);
end

denom_exp_u = denom_exp_u + 1;
prob_i = exp_u./denom_exp_u(mktid,:);

prob_i(isinf(expmu)) = 1;


meanprob = mean(prob_i, 2);
